Bright Solitary Waves in Malignant Gliomas 
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We put forward a nonlinear wave model describing the fundamental physio-pathologic features of 
an aggressive type of brain tumors: glioblastomas. Our model accounts for the invasion of normal 
tissue by a proliferating and propagating rim of active glioma cancer cells in the tumor boundary 
and the subsequent formation of a necrotic core. By resorting to numerical simulations, phase space 
analysis and exact solutions, we prove that bright solitary tumor waves develop in such systems. 
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Introduction.- Solitons arc localized wave-packets able 
to maintain their shape and speed when propagating in 
different media and under mutual collisions. The exis- 
tence of such particlc-like waves interacting elastically is 
typical of integrable nonlinear equations [l| but they also 
arise in a broader set of physical systems described by 
nonintegrable wave equations, displaying richer spatial 
interactions with their sustaining media |2j and complex 
collision scenarios Q. In the last years there has been 
an increased interest in the application of the concepts 
and tools from nonlinear physics to biology and medicine 
were nonlinear waves can potentially appear How- 
ever, the identification of nonlinear waves in oncology 
has remained very limited Q. 

Gliomas comprise a heterogeneous group of neoplasms 
that initiate in the brain or in the spine. Glioblastoma 
multiforme (GBM) is the most common and most ag- 
gressive type of glioma with poor prognosis and survival 
ranging from 12 to 15 months after diagnosis GBMs 
are composed of a mixture of poorly differentiated neo- 
plastic astrocytes. These tumors may develop (spanning 
from 1 year to more than 10 years) from lower-grade as- 
trocytomas (World Health Organization [WHO] grade II) 
or anaplastic astrocytomas (WHO grade III), but more 
frequently, they manifest de novo, presenting after a short 
clinical history, usually less than 3 months, without any 
evidence of a less malignant precursor lesion Stan- 
dard treatments of GBMs include surgery, conformal ra- 
diotherapy and drugs such as alkylating agents and an- 
tiangiogenic therapies . 

Various models have been proposed to describe specific 
aspects of GBMs many of them based on a sim- 

ple reaction-diffusion equation: the Fischer-Kolmogorov 
(FK) equation In one-dimensional scenarios the FK 
equation has solitary wave solutions of kink- type [H, EH j 
accounting for the progression of the tumor front edge, 
but in higher dimensions its analysis must resort to 
numerical methods. Essential features of high-grade 
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FIG. 1: (a) Sagittal magnetic resonance image showing a pa- 
tient with a GBM in the occipital (posterior) lobe of the brain. 
The contrast enhanced region (marked by arrows) contains 
the active tumor cells and the darker inner part of the rim 
is the necrotic tissue, (b) Schematic representation of the 
different areas. Arrows indicate the tumor and necrotic core 
growth directions. 



gliomas are neglected by the FK model, such as (i) the 
formation of a core of necrotic tissue (see Fig.[TJ), respon- 
sible for the intracranial deformation that may lead to 
death and (ii) the interaction of the tumor with adja- 
cent normal tissue. More elaborated approaches include 
additional details in (only) some of the intervening mech- 
anisms (IlT - flBI but lack sufficient biological information 
on the parameters and other processes that impair the 
model's predictive capability. 

There is thus a need for models accounting for the cru- 
cial features of GBM dynamics (see Fig.[T]) without incor- 
porating excessive details on any of the -often unknown- 
specific processes. Ideally, such models should be simple 
enough to allow for some quantitative understanding, e.g. 
using tools of nonlinear wave theory [l[ . In this letter we 
present a model capturing the key features seen in real 
GBMs and enabling us to carry out a full theoretical anal- 
ysis. The model predicts the existence of bright solitary 
waves, bright solitons, for the spreading tumor cell front. 

The model.- We will consider a simple description of 
the invasion phenomenon in high-grade gliomas by incor- 
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porating the interaction of three relevant (nonncgative) 
densities. The tumor cell density, denoted by a function 
u(x,t), the adjacent normal cells v(x, t), and a necrotic 
core density w(x,t), so that our model reads as 

du 



at 

dv 

~dt 

dw 



= DAu + p{u* 
= —J-{u,v,w). 
= F{u, v, w) + 



w)u 



au . 



(la) 
(lb) 
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where A = X)7=i ®* /® x j an d the boundary conditions 
are u(±oo) = 0, u(±oo) = v*,u;(±oo) = 0. For the tu- 
mor to grow it is necessary that pu* > a. The standard 
FK equation is recovered by setting v = w = 0, how- 
ever, in order to properly describe the observed clinical 
phenomenology, we include both the population of the 
normal tissue and the developing necrotic core, this be- 
ing a distinctive feature of GBMs (see Fig. [1]). Here, we 
will focus in the one-dimensional version of Eqs. ([TJ but 
the shown phenomena persist in higher dimensions. 

In Eqs. (p}, tumor cell spreading is incorporated using 
a standard Fickian diffusion mechanism. This is the dif- 
fusion employed in most of the continuous mathematical 
models of cell motility. Diffusion phenomena in gliomas 
should probably be governed by more complicated frac- 
tional (anomalous) diffusion [l7[ or other more elaborate 
terms [l8[ to account for the high infiltration observed 
in this type of tumors [ijj and the fact that cells do 
not behave like purely random walkers and may actu- 
ally remain immobile for a significant amount of time 
before compelled to migrate to a more favorable place. 
The nonlinear term in Eq. (|la|) corresponds to prolifera- 
tion under a competition for space between the different 
densities. A tumor cell death term is added to include 
the fact that tumor cells, although generally lacking the 
programmed cell death mechanisms [23], may succumb 
because of their competition with the immune system in 
normal tissue, hypoxia and acidosis in the high density 
tumor areas, and deficiency of nutrients and physical sup- 
port in the necrotic core. In average, the characteristic 
tumor cell life time is 1/a. 

In Eq. (fib)) we have represented the cell loss due to 
the interaction with the tumor by means of an arbitrary 
form IF(u,v,w) depending on all the densities. The de- 
tails of the interaction may be very complicated. The 
mechanisms of cell death combine the acidosis generated 
by the anomalous metabolism of the tumor cells [2lj , the 
competition for nutrients with the glioma cells, the mod- 
ification of the microenvironment [22|, and other effects. 
In agreement with physiological data, we will assume that 
the normal brain tissue does not proliferate. Finally, we 
have assumed in Eq. (fTcl) that the space occupied by the 
necrotic core is the same compartment occupied by the 
original cells and grows at the expense of the other two 
compartments. More elaborate models could include a 



reduction coefficient to account for the shrinkage of the 
cells and/or the destruction of their cytoplasm and the 
release of the cellular content to the necrotic area. 

Numerical simulations.- To elucidate the typical dy- 
namics of Eqs. (fTJ we display the results of numerical 
simulations in Figs. [2] and [3] For this particular case we 
have chosen a contact interaction term between the tu- 
mor and the normal tissue of the form J-(u, v,w) — 7 u v. 
The used parameters come from clinically observed val- 
ues for the diffusion coefficient D = 0.02 mm 2 /day, prolif- 
eration rate p = 0.5 day -1 , death rate a = 1/30 day -1 , 
and invasion parameter 7 = 0.25 day -1 . Initial data 
are taken as u(x, 0) = 0.1sech(5x), v(x, 0) = = 0.4, 
w(x,0) = (all in units of u*) corresponding to a small 
dysplasia in a bed of normal cells. 

Initially (t = 15 days in Fig. the tumor develops 
embedded in the environment of normal cells and with 
space for proliferation. As space saturates because of the 
accumulation of dead tissue and tumor cells, the cells 
can no longer survive and start depleting the tumor at 
its initial location (see e.g. t = 60 days in Fig. [2]). The 
tumor then generates a bright soliton that propagates in 
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FIG. 2: (color online) Simulations of Eqs. |T} in one spa- 
tial dimension leading to the formation of a bright soliton. 
Parameter values are: D = 0.02 mm 2 /day, p — 0.5 day -1 , 
a — 1/30 day -1 , and 7 = 0.25 day -1 . Shown are the profiles 
of the density of the different densities (a) tumor u(x,t), (b) 
normal tissue v(x,t), and (c) necrotic core w(x,t), for times 
t — 15, 60, 120 and t = 180 days. All of the densities are 
measured in units of u* . 
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FIG. 3: (color online) Simulations of Eqs. (JTJ in one spa- 
tial dimension leading to the formation of a bright soliton. 
Parameter values are as in Fig. [2] Evolution of the to- 
tal cell number (a), center (b), and speed (c) of the soli- 
ton traveling to the left, defined as N(t) = f_ u(x,t)dx, 
X(t) = f^^xu^jtydx/N^), and dX/dt, respectively. 



the normal tissue together with a kink in the normal tis- 
sue and an anti-kink in the necrotic core (c.f. Fig. [2]). 
After the initial transient the tumor mass stabilizes [see 
e.g. Fig.^a)] following the classical Gompertzian growth 
curve but now with nontrivial spatial effects incorpo- 
rated. Finally, as seen in Fig. [3{c), the speed of the 
soliton stabilizes at a constant value fundamentally dif- 
ferent from the one from the FK model and dependent 
on a [Fig. Ufa)]. As with the FK model, faster solitons 
from initial moving data may be possible, but here we 
focus our attention on the minimal speed solutions still 
arising from initial data. 

Theory.- It is convenient to introduce the new func- 
tions {/(C,t) = u/u* and S(( 7 t) = 1-/3 + [v + w] /u*, 
together with the rescaled variables C = pu*/D and 
t = pu*t. Then, Eqs. ([T]) can be cast as 



U T = U a -(U + S)U , 
S T = 0U, 



(2a) 
(2b) 



with /? = a/pu* < 1 and subscripts will henceforth de- 
note partial differentiation. Let us first notice that, sur- 
prisingly, the precise form of the tumor-host interaction 
term T is not relevant for the U-S dynamics. Instead, 
the global function S incorporating both the necrotic tis- 
sue and the normal cells accounts for the effect of the 
peritumoral environment. Secondly, because of Eq. (|2b[) . 
the non-tumoral density at any given point always in- 
creases which implies that after the tumor wave crosses 
a region, the necrotic core contains a higher density than 
the normal stroma. Thirdly, since the physically feasible 
solutions must be positive, we combine Eqs. @ into a 
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FIG. 4: (color online) Bright solitons. (a) Velocity depen- 
dence on a/pu*. The other parameters and initial data are 
as in Fig. [2] (b) Profiles from Eq. Q (solid curve) and 
the explicit solution given by Eq. ([7| (dashed curve) for 
P = a/pu, = 0.4, c = 3\/pu t D and cj>o = 0.5m,. 



single (and more general) equation for the tumor density 



= -pu. 



(3) 



Short time limit.- When < r <C 1//3, the right-hand- 
side term in Eq. §3§ can be neglected (no significant tu- 
mor cell death has occurred yet). The resulting equa- 
tion U T — U^q + U 2 = possesses self-similar solutions 
of the form U((,t) — ip (£t -1 / 2 ) /r, where ip satisfies 
2 (tp m + tp — tp 2 ) + r]ip v = and 77 = (t -1 / 2 . For U to 
remain constant in the limit t — > + , it must be the case 
that ip — > as 77 — > 00. Hence, we can disregard the non- 
linear term, and consider 2 {ip m + p) + rjip^ = 0, whose 
physically meaningful solution is tp = tp \if\ cxp(— ?/ 2 /4). 
This profile is consistent with the initial ones observed in 
the numerical simulations [sec Fig.[2Ja)]. 

Solitary waves.- Our numerical simulations show that 
two counter-propagating tumor wave fronts develop for 
i » 1/a. In this long time limit, both fronts acquire 
a characteristic traveling solitary shape (they no longer 
interact). Wc will look for such solutions of Eq. ([3]) in the 
form U(C — ct) = 4>{rj), where c is their velocity (positive 
or negative for right- or left-moving fronts, respectively). 
These solutions will be required to satisfy <fi < 1 — p, 
together with <\> —J- and <pc —> as [77! — > 00. Upon 
substitution in Eq. ([3]) and defining £ = 77/c, we arrive at 
a third order nonlinear autonomous differential equation 



i(^-#^) = 0. (4) 



The existence of traveling waves of Eq. ((4]) as a function 
of P and c is a mathematical bifurcation problem beyond 
the scope of this paper. Our numerical calculations re- 
veal that: (i) the minimum (absolute value) speed c mlIl 
for which positive solutions <f> < 1 — P exist decreases as 
P — > and can attain values c mm < 2y / pu*D. This is in 
contrast with the FK equation for which the minimum 
speed must always fulfill c mnl > 2yJ pu^D Therefore, 
our model [see Fig. EJa)] can account for slower growing 
gliomas than those predicted by the standard FK equa- 
tion [l(J. (ii) For \c\ > 5\/pu*D, the last two terms in 
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Eq. (j4| can be neglected for any 
simpler equation 



< 1 



P4? -4- + = o . 



to yield the 



(5) 



Phase portrait analysis.- Bright soliton solutions to 
Eq. ([5]) correspond to homoclinic paths in the phase 
plane. Defining %p(<p) = <f>£, Eq. {[5J transforms into the 
orbit equation = <f> + ip/tj) — (3(f) 2 /ip, whose solution, 
with x = W0; i s X + /?log \P — x\ = C + <f>. The criti- 
cal points along the orbits obey <f) = x 2 HP — x) giving 
rise to two distinct regions for <j> > 0. The region x > P 
contains no localized solutions since the orbit ip{<fi) is a 
monotonic increasing function. The region x < P does 
contain localized solutions if the integration constant C is 
bounded from above. This can be shown by considering 



f(x) = X + Plog(i3-x) 



P-x 



-c. 



(6) 



which presents a global maximum at x = 0- I n addition, 



it holds that lim^-oo /(%) = lim 



f(x) 



-CO. 



Thus, the existence of two different roots (with opposite 
signs) is guaranteed whenever C < p log/3. This proves 
the existence of fast bright soliton solutions. A complete 
theoretical analysis will be reported elsewhere. 

Explicit soliton solutions.- As a final comment, let us 
mention that the approximate explicit form of bright soli- 
ton solutions from Eq. ([5]) can be obtained by neglecting 
the third term in Eq. (|S|). This term introduces a small 
asymmetry in the profile. The found expression is 



0(0 = 0o sech 2 



(7) 



with constants 0o and £o- Fig. 0|b) compares the ex- 
act numerical profile from Eq. (01 and the one predicted 
by Eq. (JT]). The non-tumoral density S can be obtained 
by integrating Eq. (|2b|) and exhibits a kink-like shape. 
Within the necrotic compartment, this is the typical 
contrast-enhanced region seen in the magnetic resonance 
images. The tumor cell number N(t) = J a ^ oo u(x,t)dx 
for the left-moving soliton follows from Eq. ([7} 



N(t) 



a 



1 + tanh 



' 4>opu*a 



(8) 



which, for long times, gives a saturation similar to the 
one depicted in Fig. EJa). 

Conclusions.- We have developed a simple model of 
glioblastoma progression incorporating the normal tissue, 
tumor cells and the necrotic core. Our theoretical study 
displays many of the signatures of aggressive malignant 
gliomas with bright solitons acting as attractors of the 
tumor-host dynamics that can be compared with the ob- 
served phenomenology. The two and three dimensional 
versions of our starting Eqs. (TTJ can easily accommodate 



the effect of ionizing radiation both on the tumor and on 
the normal tissue [23[ . This could provide a useful tool, if 
combined with intensity-modulated radiation therapy, to 
simulate the outcome of different dose painting scenarios. 
We hope that our model will stimulate the use of tech- 
niques from nonlinear physics and nonlinear wave theory 
to further understand key aspects of tumor growth. 
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